#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Process the LCMP Data

#Download LCMAP Landcover Data
#put it into a huc 12 and year panel

#load packages
library(raster)
library(sf)
library(exactextractr)
library(tidyverse)
library(mapview)
library(data.table)

#open 2018
setwd("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")

lc2018<-raster("Data/LCMAP/LCMAP_CU_2018_V13_LCPRI.tif")

lc2018@data@values

mapview(lc2018, legend=TRUE)

#there are 8 land cover classifications
#1 developed
#2 cropland
#3 grass/shrub
#4 tree cover
#5 water
#6 wetland
#7 ice/snow
#8 barren

#legend colors
#https://imagecolorpicker.com/en 
legend<-c("developed", "cropland", "grass/shrub", "tree cover",
          "water", "wetland", "ice/snow", "barren")

colors<-c("#ff2e2e", "#bf8d5a", "#e7f0d3", "#5c906a", 
          "#0071e2", "#b4daff", "#ffffff", "#b4b0a5")


mapview(lc2018, zcol="value", at=c(1,2,3,4,5,6,7,8), 
        col.regions=colors, na.color="#ededed", legend=FALSE)


##get aggregated land use of each huc 12. 
#HUC 12
MARB<-st_read("Data/WatershedBoundaries/Processed/MARB/MARB.shp")

crs(lc2018)
crs(MARB)

MARB<-st_transform(MARB, st_crs(lc2018))
crs(MARB)

###############################################################################
#Extract

#extract land cover in 2018
landcov_huc_2018 <- exact_extract(lc2018, MARB, function(df) {
  df %%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2018)

#save it
save(landcov_huc_2018, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2018.RData")

#2017
lc2017<-raster("Data/LCMAP/Raw/LCMAP_CU_2017_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2017))

landcov_huc_2017 <- exact_extract(lc2017, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2017)

#save it
save(landcov_huc_2017, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2017.RData")

#get 2016
lc2016<-raster("Data/LCMAP/Raw/LCMAP_CU_2016_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2016))

landcov_huc_2016 <- exact_extract(lc2016, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2016)
save(landcov_huc_2016, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2016.RData")

#2015
lc2015<-raster("Water quality/Data/LCMAP/Raw/LCMAP_CU_2015_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2015))

landcov_huc_2015 <- exact_extract(lc2015, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2015)
save(landcov_huc_2015, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2015.RData")

#2014
lc2014<-raster("Data/LCMAP/Raw/LCMAP_CU_2014_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2014))

landcov_huc_2014 <- exact_extract(lc2014, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2014)
save(landcov_huc_2014, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2014.RData")

#2013
lc2013<-raster("Data/LCMAP/Raw/LCMAP_CU_2013_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2013))

landcov_huc_2013 <- exact_extract(lc2013, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2013)
save(landcov_huc_2013, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2013.RData")

#2012
lc2012<-raster("Data/LCMAP/Raw/LCMAP_CU_2012_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2012))

landcov_huc_2012 <- exact_extract(lc2012, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2012)
save(landcov_huc_2012, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2012.RData")

rm(landcov_huc_2012)
rm(lc2012)

#2011
lc2011<-raster("Data/LCMAP/Raw/LCMAP_CU_2011_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2011))

landcov_huc_2011 <- exact_extract(lc2011, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2011)
save(landcov_huc_2011, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2011.RData")

rm(lc2011)
rm(landcov_huc_2011)

#2010
lc2010<-raster("Data/LCMAP/Raw/LCMAP_CU_2010_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2010))

landcov_huc_2010 <- exact_extract(lc2010, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2010)
save(landcov_huc_2010, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2010.RData")

#2009
lc2009<-raster("Data/LCMAP/Raw/LCMAP_CU_2009_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2009))

landcov_huc_2009 <- exact_extract(lc2009, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2009)
save(landcov_huc_2009, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2009.RData")

rm(landcov_huc_2009)
rm(lc2009)

#2008
lc2008<-raster("Data/LCMAP/Raw/LCMAP_CU_2008_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2008))

landcov_huc_2008 <- exact_extract(lc2008, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2008)
save(landcov_huc_2008, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2008.RData")

rm(landcov_huc_2008)
rm(lc2008)

#2007
lc2007<-raster("Data/LCMAP/Raw/LCMAP_CU_2007_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2007))

landcov_huc_2007 <- exact_extract(lc2007, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2007)
save(landcov_huc_2007, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2007.RData")

rm(landcov_huc_2007)
rm(lc2007)

#2006
lc2006<-raster("Data/LCMAP/Raw/LCMAP_CU_2006_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2006))

landcov_huc_2006 <- exact_extract(lc2006, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2006)
save(landcov_huc_2006, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2006.RData")

rm(landcov_huc_2006)
rm(lc2006)

#2005
lc2005<-raster("Data/LCMAP/Raw/LCMAP_CU_2005_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2005))

landcov_huc_2005 <- exact_extract(lc2005, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2005)
save(landcov_huc_2005, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2005.RData")

rm(landcov_huc_2005)
rm(lc2005)

#2004
lc2004<-raster("Data/LCMAP/Raw/LCMAP_CU_2004_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2004))

landcov_huc_2004 <- exact_extract(lc2004, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2004)
save(landcov_huc_2004, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2004.RData")

rm(landcov_huc_2004)
rm(lc2004)

#2003
lc2003<-raster("Data/LCMAP/Raw/LCMAP_CU_2003_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2003))

landcov_huc_2003 <- exact_extract(lc2003, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2003)
save(landcov_huc_2003, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2003.RData")

rm(landcov_huc_2003)
rm(lc2003)

#2002
lc2002<-raster("Data/LCMAP/Raw/LCMAP_CU_2002_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2002))

landcov_huc_2002 <- exact_extract(lc2002, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2002)
save(landcov_huc_2002, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2002.RData")

rm(landcov_huc_2002)
rm(lc2002)

#2001
lc2001<-raster("Data/LCMAP/Raw/LCMAP_CU_2001_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2001))

landcov_huc_2001 <- exact_extract(lc2001, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2001)
save(landcov_huc_2001, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2001.RData")

rm(landcov_huc_2001)
rm(lc2001)

#2000
lc2000<-raster("Data/LCMAP/Raw/LCMAP_CU_2000_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc2000))

landcov_huc_2000 <- exact_extract(lc2000, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_2000)
save(landcov_huc_2000, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2000.RData")

rm(landcov_huc_2000)
rm(lc2000)

#1999
lc1999<-raster("Data/LCMAP/Raw/LCMAP_CU_1999_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1999))

landcov_huc_1999 <- exact_extract(lc1999, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1999)
save(landcov_huc_1999, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1999.RData")

rm(landcov_huc_1999)
rm(lc1999)

#1998
lc1998<-raster("Data/LCMAP/Raw/LCMAP_CU_1998_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1998))

landcov_huc_1998 <- exact_extract(lc1998, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1998)
save(landcov_huc_1998, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1998.RData")

rm(landcov_huc_1998)
rm(lc1998)

#1997
lc1997<-raster("Data/LCMAP/Raw/LCMAP_CU_1997_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1997))

landcov_huc_1997 <- exact_extract(lc1997, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1997)
save(landcov_huc_1997, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1997.RData")

rm(landcov_huc_1997)
rm(lc1997)

#1996
lc1996<-raster("Data/LCMAP/Raw/LCMAP_CU_1996_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1996))

landcov_huc_1996 <- exact_extract(lc1996, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1996)
save(landcov_huc_1996, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1996.RData")

rm(landcov_huc_1996)
rm(lc1996)

#1995
lc1995<-raster("Data/LCMAP/Raw/LCMAP_CU_1995_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1995))

landcov_huc_1995 <- exact_extract(lc1995, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1995)
save(landcov_huc_1995, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1995.RData")

rm(landcov_huc_1995)
rm(lc1995)

#1994
lc1994<-raster("Data/LCMAP/Raw/LCMAP_CU_1994_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1994))

landcov_huc_1994 <- exact_extract(lc1994, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1994)
save(landcov_huc_1994, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1994.RData")

rm(landcov_huc_1994)
rm(lc1994)

#1993
lc1993<-raster("Data/LCMAP/Raw/LCMAP_CU_1993_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1993))

landcov_huc_1993 <- exact_extract(lc1993, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1993)
save(landcov_huc_1993, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1993.RData")

rm(landcov_huc_1993)
rm(lc1993)

#1992
lc1992<-raster("Data/LCMAP/Raw/LCMAP_CU_1992_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1992))

landcov_huc_1992 <- exact_extract(lc1992, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1992)
save(landcov_huc_1992, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1992.RData")

rm(landcov_huc_1992)
rm(lc1992)

#1991
lc1991<-raster("Data/LCMAP/Raw/LCMAP_CU_1991_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1991))

landcov_huc_1991 <- exact_extract(lc1991, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1991)
save(landcov_huc_1991, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1991.RData")

rm(landcov_huc_1991)
rm(lc1991)

#1990
lc1990<-raster("Data/LCMAP/Raw/LCMAP_CU_1990_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1990))

landcov_huc_1990 <- exact_extract(lc1990, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1990)
save(landcov_huc_1990, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1990.RData")

rm(landcov_huc_1990)
rm(lc1990)

#1989
lc1989<-raster("Data/LCMAP/Raw/LCMAP_CU_1989_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1989))

landcov_huc_1989 <- exact_extract(lc1989, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1989)
save(landcov_huc_1989, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1989.RData")

rm(landcov_huc_1989)
rm(lc1989)

#1988
lc1988<-raster("Data/LCMAP/Raw/LCMAP_CU_1988_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1988))

landcov_huc_1988 <- exact_extract(lc1988, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1988)
save(landcov_huc_1988, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1988.RData")

rm(landcov_huc_1988)
rm(lc1988)

#1987
lc1987<-raster("Data/LCMAP/Raw/LCMAP_CU_1987_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1987))

landcov_huc_1987 <- exact_extract(lc1987, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1987)
save(landcov_huc_1987, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1987.RData")

rm(landcov_huc_1987)
rm(lc1987)

#1986
lc1986<-raster("Data/LCMAP/Raw/LCMAP_CU_1986_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1986))

landcov_huc_1986 <- exact_extract(lc1986, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1986)
save(landcov_huc_1986, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1986.RData")

rm(landcov_huc_1986)
rm(lc1986)

#1985
lc1985<-raster("Data/LCMAP/Raw/LCMAP_CU_1985_V13_LCPRI.tif")
MARB<-st_transform(MARB, st_crs(lc1985))

landcov_huc_1985 <- exact_extract(lc1985, MARB, function(df) {
  df %>%
    mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
    dplyr::group_by(ID, value) %>%
    dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = c('ID', 'huc12'), progress = TRUE)

list(landcov_huc_1985)
save(landcov_huc_1985, file="Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1985.RData")

rm(landcov_huc_1985)
rm(lc1985)


################################################################################
#Putting it all together.
#add the years and bind the rows

#open files
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1985.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1986.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1987.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1988.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1989.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1990.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1991.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1992.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1993.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1994.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1995.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1996.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1997.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1998.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_1999.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2000.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2001.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2002.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2003.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2004.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2005.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2006.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2007.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2008.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2009.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2010.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2011.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2012.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2013.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2014.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2015.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2016.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2017.RData")
 load("Data/LCMAP/Processed/yearly_huc12/LCMAP_Landcov_huc_2018.RData")

#add year variable
landcov_huc_1985$year<-1985
landcov_huc_1986$year<-1986
landcov_huc_1987$year<-1987
landcov_huc_1988$year<-1988
landcov_huc_1989$year<-1989
landcov_huc_1990$year<-1990
landcov_huc_1991$year<-1991
landcov_huc_1992$year<-1992
landcov_huc_1993$year<-1993
landcov_huc_1994$year<-1994
landcov_huc_1995$year<-1995
landcov_huc_1996$year<-1996
landcov_huc_1997$year<-1997
landcov_huc_1998$year<-1998
landcov_huc_1999$year<-1999
landcov_huc_2000$year<-2000
landcov_huc_2001$year<-2001
landcov_huc_2002$year<-2002
landcov_huc_2003$year<-2003
landcov_huc_2004$year<-2004
landcov_huc_2005$year<-2005
landcov_huc_2006$year<-2006
landcov_huc_2007$year<-2007
landcov_huc_2008$year<-2008
landcov_huc_2009$year<-2009
landcov_huc_2010$year<-2010
landcov_huc_2011$year<-2011
landcov_huc_2012$year<-2012
landcov_huc_2013$year<-2013
landcov_huc_2014$year<-2014
landcov_huc_2015$year<-2015
landcov_huc_2016$year<-2016
landcov_huc_2017$year<-2017
landcov_huc_2018$year<-2018

#bind data
landcov_huc<-dplyr::bind_rows(landcov_huc_1985,
                              landcov_huc_1986,
                              landcov_huc_1987,
                              landcov_huc_1988,
                              landcov_huc_1989,
                              landcov_huc_1990,
                              landcov_huc_1991,
                              landcov_huc_1992,
                              landcov_huc_1993,
                              landcov_huc_1994,
                              landcov_huc_1995,
                              landcov_huc_1996,
                              landcov_huc_1997,
                              landcov_huc_1998,
                              landcov_huc_1999,
                              landcov_huc_2000,
                              landcov_huc_2001,
                              landcov_huc_2002,
                              landcov_huc_2003,
                              landcov_huc_2004,
                              landcov_huc_2005,
                              landcov_huc_2006,
                              landcov_huc_2007,
                              landcov_huc_2008,
                              landcov_huc_2009,
                              landcov_huc_2010,
                              landcov_huc_2011,
                              landcov_huc_2012,
                              landcov_huc_2013,
                              landcov_huc_2014,
                              landcov_huc_2015,
                              landcov_huc_2016,
                              landcov_huc_2017,
                              landcov_huc_2018)

rm(landcov_huc_1985,
   landcov_huc_1986,
   landcov_huc_1987,
   landcov_huc_1988,
   landcov_huc_1989,
   landcov_huc_1990,
   landcov_huc_1991,
   landcov_huc_1992,
   landcov_huc_1993,
   landcov_huc_1994,
   landcov_huc_1995,
   landcov_huc_1996,
   landcov_huc_1997,
   landcov_huc_1998,
   landcov_huc_1999,
   landcov_huc_2000,
   landcov_huc_2001,
   landcov_huc_2002,
   landcov_huc_2003,
   landcov_huc_2004,
   landcov_huc_2005,
   landcov_huc_2006,
   landcov_huc_2007,
   landcov_huc_2008,
   landcov_huc_2009,
   landcov_huc_2010,
   landcov_huc_2011,
   landcov_huc_2012,
   landcov_huc_2013,
   landcov_huc_2014,
   landcov_huc_2015,
   landcov_huc_2016,
   landcov_huc_2017,
   landcov_huc_2018)

value<-c(1:8)
key<-as.data.frame(value)
key$landcover<-c("developed", "cropland", "grass_shrub", "tree_cover",
          "water", "wetland", "ice_snow", "barren")

landcov_huc<-dplyr::left_join(landcov_huc, key, by="value")
landcov_huc<- landcov_huc %>% relocate(landcover, .after=value)

#ID and huc12 key
huc12key<-subset(as.data.frame(MARB), select=c(ID, huc12))
rm(MARB)

landcov_huc<-dplyr::left_join(landcov_huc, huc12key, by="ID")
landcov_huc<- landcov_huc %>% relocate(huc12, .after=ID)

summary(landcov_huc)

landcov_huc<-subset(landcov_huc, select=c(year, huc12, landcover, freq))
rm(key)
rm(huc12key)

#turn into wide pane#turn from long to wide
landcov_agg_wide<-dcast(landcov_huc, huc12 + year ~ landcover, value.var="freq")
summary(landcov_agg_wide)

#if its NA, replace with 0. 
landcov_agg_wide<-landcov_agg_wide %>% replace(is.na(.), 0)
summary(landcov_agg_wide)

lcmap_huc12<-landcov_agg_wide[,c(1:10)]

table(lcmap_huc12$year)
summary(lcmap_huc12)

#save the data
save(lcmap_huc12, file="Data/LCMAP/Processed/lcmap_landcover_huc12_panel.RData")
